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Abstract 

We investigate the influence of boundary slip velocity in Newtonian 
fluids at finite Reynolds numbers. Numerical simulations with Lattice 
Boltzmann method (LBM) and Finite Differences method (EDM) are 
performed to quantify the effect of heterogeneous boundary conditions 
on the integral and local properties of the flow. Non linear effects 
are induced by the non homogeneity of the boundary condition and 
change the symmetry properties of the flow inducing an overall mean 
flow reduction. To explain the observed drag modification, recipro- 
cal relations for stationary ensembles are used, predicting a reduction 
of the mean flow rate from the creeping flow to be proportional to 
the fourth power of the friction Reynolds number. Both numerical 
schemes are then validated within the theoretical predictions and re- 
veal a pronounced numerical efficiency of the LBM with respect to 
EDM. 

1 Introduction 

The growing interest in fluid properties and mass transfer at the micro and 
nanoscale [H |2] has recently produced new research themes and questions. 
Above all, the clear understanding of interfacial phenomena and wetting 
properties jHl Ej , roughness effects [3 El Ej and surface nanobubbles jH] (see 
also for an exhaustive review) is constantly providing new perspectives on 



1 



how to interpret boundary conditions for fluid flows confined to the micro- 
scale and even below. In fact, the failure of the classical no-slip boundary 
condition in hydrodynamics [lOj, is now predicted by a series of experiments 
[TTl 1121 CSl mi and numerical studies [T31 HE] ■ An appealing explanation 
for this observed slippage is the formation of gas pockets (bubbles) between 
the liquid and the solid ^7j, in such a way that they can provide a zero shear 
stress condition for the flow and modify considerably its friction properties. 
One would then be interested to understand which are the microscopic pa- 
rameters at the onset of this slip motion and that, more generally, determine 
the boundary stability of these bubbles. Such an ab-initio study should cor- 
rectly be addressed in Molecular Dynamics simulations (MD) fSJ |B] where 
the problem is attacked from the atomistic point of view, integrating New- 
ton's equation for a set of molecules interacting through Lennard- Jones po- 
tentials. The main drawback of this approach is (by definition) its inability 
to describe spatial fluctuations larger than the intermolecular potential and 
temporal scales larger than a few milliseconds. This obviously prevents a 
clear understanding of the slip effects on the flow properties on larger scales, 
an issue that must be addressed through continuum JHl UH 1201 121] descrip- 
tion based on the Navier-Stokes equations with slip boundary conditions or, 
alternatively, with recently used mesoscopic methods based on the Boltzmann 
equation (the Lattice Boltzmann method (LBM)) with kinetic boundary con- 
ditions [221 123 1211 • More precisely, in order to study the correct continuum 
momentum balance one should work on the continuity and Navier-Stokes 
equations: 

V • n = 

Vp 

dtu + V ■ {uu) = h uAu 

P 

where u is the velocity vector, p the internal pressure, p the fluid density 
and u the kinematic viscosity. These equations should then be provided 
with some ad hoc boundary conditions with the role of renormalizing the 
microscopic details of the fluid interactions at the wall fOl 10] • When such a 
continuum description is involved, the first non trivial control parameter on 
the flow is the Reynolds Number, 

V 

expressing the ratio between the importance of advective terms with respect 
to viscous ones, being L,!!,^ a. typical macroscopic length, macroscopic ve- 



2 



locity and fluid viscosity respectively. Small scale systems imply small Re 
{Re -C 1), the consequence being that non linear terms are considered neg- 
ligible and consequently ruled out. This means that when the boundary 
conditions are declared, we have to face a linear problem that in some cases 
is feasible from the mathematical point of view. In particular, the effective 
slip properties of pressure driven Stokes flows over heterogeneous surfaces 
made up of mixed no slip and free shear walls (hereafter NSW and FSW) 
have been studied 

1„ 1 fdP\ 
uAu = - \/p H — — 

p p \dx J 

■U|| = (NSW) 
dnU\\ = (FSW) 

being ^ an external pressure pumping (required to obtain a nonzero mass 
flow rate) in the direction e^,., U\\ the tangential component of the fluid ve- 
locity at the wall and n the wall normal directed outward in the bulk fluid. 
The physical idea is to translate the presence of gas pockets in a suitable 
boundary condition for the hydrodynamic flelds and then study the drag 
modiflcations in the fluid as a function of the degrees of freedom describing 
the surface heterogeneities. This problem has been addressed for the flrst 
time by Philip nearly 30 years ago |19| 1201 for the case of longitudinal strips 
of free shear (with respect to the flow direction). Using conformal mapping 
Philip related the mass flow rate gain due to the slip motion to the geometry 
of the boundary condition. Only more recently Stone and Lauga JHl have 
proposed a similar approach for pressure driven Stokes flows with transversal 
strips of free shear. In these models, non linear terms are supposed very small 
(formally Re —>■ 0) but, due to the surface heterogeneities, it is not clear if 
and how they influence the fluid by changing Re. In fact, due to the non 
homogeneity in the boundary condition, non subtle effects may be induced 
by the boundary and make the non linear terms play a role, an effect that 
would have no counterpart in laminar homogeneous flows. To clarify this 
point, in this paper we will study the overall sensitivity of the fluid with re- 
spect to these effects. We will carry out a complete study of pressure driven 
laminar flows with mixed boundary conditions of free shear and no slip for 
finite Reynolds numbers : 

V ■ [UUj = U/AU Vp — — 

p p \dx J 
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nil = (NSW) 

dnU\\ = (FRW). 

Numerical Simulations with Finite Difference methods (FDM) and Lattice 
Boltzmann methods (LBM) are used firstly to check the analytical estimates 
proposed by Stone and Lauga (for Re = 0) in the case of transversal strips 
and then to characterize the flow properties at Re > 0. Interesting flow 
behaviors are observed for finite Re and a non trivial interplay between the 
boundary condition and non linear terms is considered both theoretically and 
numerically. In particular, the mass flow rate, (tt), for finite Reynolds will 
differ from its creeping flow [Re = 0) counterpart, (u), in an exact way: 

{u)-{u) = {iuu):iVu)). 

Employing finite Reynolds perturbation theory, for small Reynolds numbers, 
we will show that the observed drag modification, even if triggered by the 
boundary, has a scaling law behavior with respect to the bulk Reynolds 
number 

{u)/{u) - 1 ~ i?el 

In order to validate this scaling law behavior and its range of applicability, all 
the observed effects are quantitatively studied with numerical simulations. 
The result is that up to i?e ~ 1 the predicted scaling law behavior is satisfied 
with high accuracy in both numerical schemes and the friction properties of 
the system are influenced by the presence of non linear terms in such a way 
that the flow differs from its creeping flow counterpart both qualitatively 
and quantitatively. In fact soon after non linear effects are induced by the 
boundary conditions, the flow is still laminar but its local symmetry proper- 
ties change as an overall drag enhancement is induced. 
The paper is organized as follows: in section |21 the physical problem is math- 
ematically formulated and we will give the correct background to apply flnite 
Reynolds perturbation theory that is the subject of section|21 In sectionlHwe 
give a brief review of the numerical procedures used. All numerical results 
are discussed in section |5] and conclusions will follow in section IHl 

2 Formulation of the problem 

In this section we formulate the problem under consideration from the math- 
ematical point of view. We will refer to a wall-bounded flow where x, y and 
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z denote respectively the streamwise, spanwise and wall-normal coordinates. 
The governing equations are expressed in the dimensional form: 

V* • w* = 0, (1) 



1 / dP*\ 

V* • = V* • a* + - (^-^ j e., (2) 

where u*{= {u*,v*,w*)) is the velocity vector, cr* the stress tensor, p* the 
density of fluid, —dP*/dx* the driving pressure gradient, and the unit 
vector in the ,t— direction. In the above equations, the superscript * has 
been introduced in order to indicate dimensional quantities and to distinguish 
them from dimensionless quantities (without superscript) to be introduced 
later. 

For the incompressible Newtonian fluid, the stress tensor is written as 

a* ^ -—I + 2u*S*, (3) 

where p* is the deviation of the pressure from the driving component, / the 
unit dyadic, u* the kinematic viscosity and S* the strain rate given by 

S* = l{V*u* + (V*u*f}. (4) 

Boundary conditions on the no-slip wall (NSW) and the free slip wall (FSW) 
are written as 

It* = on NSW, (5) 



(n ■ o-*) X n = on FSW, (6) 

where n is the wall-normal unit vector. Two walls are then considered at 
z — ^ and z = Lz while in the streamwise and spanwise directions, periodic 
boundary conditions are imposed with the period length of and Ly. Note 
that the relation of n • cr* • u* = is satisfied on both the no shp and free 
shear walls. In the present study, we use symmetric boundary conditions in 
the upper and lower walls and we will refer to the slip percentage, ^, as the 
ratio between the free shear area with respect to the total one. 
To proceed further, we should introduce dimensionless equations to which we 
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will apply finite Reynolds perturbation theory. To do so, let us now consider 
the friction velocity u-,- and the friction Reynolds number Re-,- defined as 



Ur = ^/T^, (7) 



fie, = MM?), (8) 

where Ty, is the skin friction averaged over the wall surface ai z — (or 
z — Lz). Since symmetric boundary conditions are imposed on z = and 
z = Lz, we obtain the relation between the skin friction and the driving 
pressure gradient: 

du\ / du\ Lz f dP\ 



(9) 



that immediately imply: 

Re oc Rel. (10) 

Now, we choose u^- and Lz/2 as the velocity and length scales for normaliza- 
tion: 

U — {u, V, w) — {u*/Ur, V*/Ur, W* / Ur) (H) 

X = (x, y, z) = (2x7L„ 2y*/Lz, 2z* jU). (12) 

This means that the two walls are located at z = and z = 2 and that the ve- 
locity components and the pressure fields satisfy the following dimensionless 
equations: 

V-w = 0, (13) 

V • (ww) = VfT (14) 
with the boundary conditions 

w = on NSW, (15) 

(n • 0-) X n = on FSW, (16) 
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where the pressure has been normahzed with the term {p*u^) and the dimen- 
sionless stress tensor is given by 

CT = -pI+ -^S. 

We also introduce the velocity vector ii and the pressure p satisfying the 
following the equations for the creeping flow with no advection: 

V ■ ii = 0, (17) 

= V ■ a- + e,, (18) 

with the boundary conditions 

w = on NSW, (19) 

(n ■ ct) X n = on FSW. (20) 

In this case, u and p represent the zeroth order approximation for the velocity 
field and stress tensor and they will be used to quantify the effects of a finite 
Reynolds number. Our starting point is the application of reciprocal relations 
|2S1 to the fields u and u and their stresses. In fact it can be shown that 

(Vw) : a- = (Vw) : a- (21) 

is an exact relation between the above fields satisfied in every point of the 
space. Taking the volume integral of the Ihs of ()21|) over the whole fluid 
region, we obtain 



d^x iVu) : cr = - d^x u ■ (V ■ cr) + 

fluid ifiuid 



/ d^x V ■ (w • cr) 

ifluid 



d'^x u ■ + (j) d^x n ■ a ■ u (22) 
fluid /wall ' 

where (■) represents the volume average operator. In the similar way, for the 
rhs of (PT|) . we obtain 

d^x (Vw) : cr = d'^x u ■ — / dc'x u ■ {{u ■ V)n} + ® d^x n ■ cr ■ u 

fluid ^ fluid ^ fluid ^wall 

-LxLzLz{u) — f d?x u ■ {{u ■ 'V)u} . 
/fluid 

(23) 
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where the second term in the rhs is 

d^a; u ■ {{u ■ V)u} = / d^x V • {u{u ■ u)} + / d^a; (uu) : (Vn) 
fluid Jfluid Jfluid 

= (b d^x n ■ {u{u ■ u)} +LxLyLz{{uu) : (Vw)). 
/wall ' V ' 

=0 

(24) 

Using (j22I) , (I2S1) and ()24|1 we then obtain the contribution of the advection to 
the change of the flow rate 

{u)-{u) = {iuu):iVu)). (25) 

Note that this expression is consistent with the identity of the Reynolds shear 
stress contribution to the friction coefficient in a turbulent flow bounded with 
no-slip walls derived in .26^ where the authors derived the following relation 



12 6{{1 - z){-u'w')) 

^'-We^ W ' ^ ^ 

where is the friction coefficient, z the wall normal coordinate normalized 



by the channel half width, Lz/2, and —u'w' the Reynolds shear stress. This 
identity indicates that the effect of the Reynolds shear stress on the skin 
friction is weighted by the factor oil — z. Substituting Vw = e^EzReril — z) 
with the no-slip walls, which contains the same weighted factor of 1 — into 
fl25|l and considering (u) = Rer/3, one can obtain the same relation. 



3 Finite Rcr Perturbation theory 

Using the perturbation method with respect to the friction Reynolds number 
Re-r under the assumption of < Rcr -C 1, the effect of the advection on 
the change of the flow rate (u) — (u) is now examined. More precisely, we 
do not directly solve the perturbation equation but we want to consider the 
following relation 

{u) - {u) = TRe"^ 

(here F is a coefficient) for small Reynolds numbers and concentrate on the 
exponent n. In our present approach, u and p in (fTSjl and ()14|) are regularly 
expended with respect to Rcr in the following form 

u = Re\u^^^ + Relu^^^ + Re^u^'^^ + (27) 
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p = Relp^"^^ + i?e^p(i) + i?e^p(2) ^ (28) 

where the properties with the overscript (0), indicate the creeping flow solu- 
tion, 

u = RerU^'^\ p = p^^\ (29) 

Note that the even order (i.e., 0{Re^)) components of the velocity vector 
are identically zero because there is not any advection or driving forcing 
to balance with the divergence of the stress tensor containing the 0{Re^) 
velocity. A detailed calculation is shown in Appendix with the final result 
that the change in the fiow rate is such that 

(n) - {u) = TRel, 

with the prefactor 

r = (ii«-{(iiW-V)n(°)}> 

that is dependent upon the perturbed velocity component but independent 
of the Reynolds number. Moreover, since (u) = Rer{u^^^) is proportional to 
Rct, the relative change of the fiow rate is such that 

{u)/{u) - 1 ~ i?e^. (30) 

Alternatively if one is interested in a scaling relation with respect to the bulk 
Reynolds number, it follows immediately from (jiup that 

{u)/{u) - 1 ~ Re\ (31) 

Note that this effect is physically induced by the boundary condition (F > 0) 
but soon after Rcr > it develops with very well defined scaling laws with 
respect to the Reynolds number. Note also that for any finite Reynolds num- 
ber, we should expect the velocity to become uniformly lower than that in 
the creeping fiow. In fact, according to Lamb |2ZI; the point in which the 
velocity field takes the minimum of the total dissipation rate is the solution 
of the creeping fiow with the prescribed boundary condition. If the fiow rate 
is fixed, the total dissipation rate is proportional to the driving pressure gra- 
dient, which equals the skin friction divided by the channel half width. This 
means that the deviation of the velocity from the creeping fiow due to the 
advective motion of fiuid enhances the total skin friction. It can be thus in- 
terpreted that a larger driving force than that in the creeping fiow is required 
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at the finite Reynolds number in order to maintain a given flow rate. Equiva- 
lently, with flxed driving pressure gradient and viscosity, the flow rate should 
reduce due to the advective motion. Therefore, (u) — {u) is negative at the 
non-zero Reynolds number, that is equivalent to (u*-^-* ■ |(n'^'^) ■ V)'U^°''}) < 
for Rcr > 0. 

4 Numerical Procedures 

In order to validate the previous theoretical analysis we will use standard 
Finite Difference methods (FDM) and Lattice Boltzmann methods (LBM). 
In the Finite Difference Method, the equations are discretized in an Eulerian 
framework on the staggered grid |2B]- The second-order scheme, i.e., the 
Adams-Bashforth method for the advection term and the Crank-Nicolson 
one for the viscous term, is used to integrate the equations in time [21] • The 
pressure is treated implicitly. The space derivatives are approximated by the 
fourth-order central difference scheme. In particular for the advection term, 
we employ the scheme by Kajishima et al. in [HO] . 

In solving the Poisson equation, we use the Fast Fourier Transform (FFT) for 
the high speed and accuracy. The three-dimensionally discretized equations 
are then reduced into the one-dimensional problem by taking the FFT in 
the homogeneous (streamwise and spanwise) directions. The reduced-order 
equation written in the heptadiagonal matrix form is directly solved. Since 
we use the staggered grid system, the velocity parallel to the wall u\\ is not 
located on the wall. In order to approximate the boundary conditions on 
the wall, the velocity at the virtual point outside the flow region is adjusted 
using the third-order Lagrange extrapolation. 

The other numerical technique used here is a mesoscopic approach based on 
the Boltzmann equations and known as Lattice Boltzmann Method (LBM) 
pn [32]. This method is a kinetic approach to fluid flows that starts from 
the Boltzmann equation in discrete form 

U{x + Atc„ t + At) - U{x, t) = - — UXx, t) - /(^-^^ {x, t)) + F,. (32) 

r 

The left hand side represents a streaming term of a probability density func- 
tion {fi{x,t)) to flnd in a given position and time {x,t) a kinetic parti- 
cle whose velocity is Cj. Here the set of velocities is properly discretized 
(cj, 2 = 0, A^) dependently on the symmetry properties of the lattice. The 
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right hand side represents a simphfied version of the standard colhsion term 
in the real Boltzmann equation: it expresses a relaxation (with character- 
istic time r) towards a local equilibrium {f^'^'^^), the equivalent of the local 
Maxwellian equilibrium in kinetic theory Finally the term Fi is an ex- 
ternal volume forcing, used to produce a constant pressure drop (dP/dx) 
in the streamwise direction. Starting from the kinetic equations (j32j) and 
coarse-graining (in the velocity space) the kinetic distributions, one obtains 
respectively a macroscopic density and velocity field as follows: 

p{x,t)=J2M^^t) (33) 

i 

uix,t) = -^J2''if^(''^^)- (34) 

Typically, to connect the above macroscopic fields to the continuum descrip- 
tion of Navier-Stokes equations one needs a separation of scale parameter, 
usually the Knudsen Number 

Kn = X/L (35) 

being A ~ \cst\ a length scale of the order of the mean free path of the kinetic 
particles (here Cs denotes the sound speed velocity), and L a macroscopic 
length. In the limit Kn -C 1 it can be shown [221 IHl] that the macroscopic 
fields fjH!?|l and (p?^ satisfy the continuum equations of fluid mechanics: 

dtp + V- (pu) = (36) 



p[dtu + {u ■ V)u] = - Vp + puAu - —e^ (37) 

being p the internal pressure and u the fluid viscosity, directly related to the 
relaxation time of the Boltzmann equation, u = c^(r — At/2) (see also 
for an exhaustive review of the method). 

The subject of boundary condition for the Lattice Boltzmann Method has 
been extensively debated in the past years. An adequate set of kinetic bound- 
ary conditions has been proposed in [21] through a lattice transcription of 
the accommodation coefficients used in continuum kinetic theory of rarefied 
gases [33] and later similar kinetic boundary condition approaches have been 
proposed on more empirical grounds [^niiSSliSZj- In particular in [22| it has 
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been shown that the introduction of a slip function at the kinetic level is 
enough to produce, for small Knudsen numbers, local macroscopic boundary 
conditions of the form: 

u\\ix) ~ Kn^^^^dnU\\{x). (38) 

The slip function s{x) represents a local degree of slippage in x, tunable 
from pure no-slip properties (s = 0) up to free shear condition (s = 1). This 
local and tunable slip properties has already been used in to match the 
macroscopic calculation proposed by Stone and Lauga in jlH] and Philip in 
[ini EOj for boundary condition made up of strips of free shear and no-slip. 

5 Numerical Results 

For the case of numerical results we concentrate on the case of transversal 
slip strips. To do so we will carry out numerical simulations with a 2d chan- 
nel with dimensions (stream-wise) and Ly (wall normal). The boundary 
condition of free shear is concentrated in a segment H on both top and lower 
walls with a slip percentage ^ immediately given by ^ = H/L^ (see figure [l}. 
The computations are made until the flow reaches a steady state. Monitoring 
the change of the kinetic energy budget and the time derivative terms is used 
to judge the steady state itself. 

Figure El shows the time evolution of the averaged kinetic energy {ul + Uy) 
as a function of time for a friction Reynolds number i?e^ = 2.245, grid mesh 
Lj. = 64, Lj^ = 84 and constant density p = 1. LBM and FDM are used 
and in both cases with the same viscosity {u = 1/6), driving pressure gra- 
dient {—dP/dx = 1.88 X 10~^) and slip percentage ^ = 0.5 . The same 
kind of plot is shown in the inset of figure El for a higher friction Reynolds 
number Rr = 10.04 and parameters = Q4:,Ly = 168, = 1/6,— dP/dx = 
4.72 X 10^^,^ = 0.5. Let us notice that the agreement on the kinetic en- 
ergy evolution is excellent with both schemes and it converges more rapidly 
towards a stationary state in the case of LBM due to the fact that in the 
kinetic scheme the Poisson equation for the pressure has not to be solved as 
in the case of FDM. In particular for the case with Rcr = 2.245 the LBM is 
5 times faster and in the case with i?e^ = 10.04 is up to 15 times faster. 
In Figure |21 we show the velocity profiles along the channels for the case 
with Rct- = 2.245 considered in figure El In the inlet region the profile of 
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Poiseuille type (parabolic profile with zero slip length) while in the middle 
of the slip region, as one can see, the local stress at the boundaries is zero 
and consequently the local viscous stress is minimized due to this free shear 
condition. As we can see not only global local properties (as it is the case of 
figure 121) but also local properties from both numerical schemes agree very 
well. 

In order to quantify the slip effects in the numerical simulations, we use the 
slip length Ls as the distance from the wall at which the extrapolated velocity 
profile is zero. To do so, we use the volume averaged streamwise velocity 



and extrapolate at zero the boundary velocity. In figure |3] the volume aver- 
aged slip length is shown as a function of the slip localization ^ for a friction 
Reynolds number Re = 2.245. The slip length has been obtained from steady 
state configurations with FDM and LBM and has been extrapolated using a 
first order accuracy and a fourth order accuracy. From the numerical results, 
the agreement with the analytical estimate and the numerical results is quite 
satisfactory: small discrepancy between the two schemes is observed for ^ 
approaching 1 but the overall trend is very well reproduced by both. 
Figure El shows the relative departure from the creeping fiow solution 1 — 
(u) I (u) for various friction Reynolds numbers. For small Reynolds number 
(u) I (u) is almost 1 while, in increasing the Reynolds number it decreases 
due to the advective motion. Within the range of friction Reynolds numbers 
shown, no turbulent eddy is observed. If the whole wall were an homogeneous 
no slip, {u)/{u) would be independent of i?e^. On the other hand, when the 
mixed boundary condition (as it is this case) is switched on, (u) / {u) is de- 
pendent on Rcr, even in the laminar regimes under consideration. According 
to our analytical investigation, for a given boundary condition (fixing ^), the 
drag modification is proportional to the fourth power of the friction Reynolds 
number and this can be verified explicitly from this plot. More precisely, re- 
lation would imply 



and to verify the correctness of this relation, together with 1 — (u) / (u) we 
plot also the right hand side of ()39|) and we observe an excellent agreement 
between the two contributions. 




1 - {u)/{u) = {{uu) : {Vu))/{u) 



(39) 
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In figure [Ul for the same set of parameters as figure we show the relative 
change of the slip length with respect to its creeping flow solution. As ex- 
pected, the same scaling behavior is observed and finite Reynolds effects 
become of the order of 10 percent soon after the friction Reynolds number is 
of the order Rcr ~ 10. 

In figures [7| and |H1 the streamwise velocity is show for two Reynolds friction 
numbers Rcr = 2.245 and Rcr = 10.04 (the same as figure Q]) . While for 
the smallest Rcr = 2.245 the profile is almost symmetric around the slip 
area, soon after Rcr is increased an asymmetry, due to the advective flow 
is observed (see figure |H1). This asymmetry is responsible for the drag mod- 
ification observed in the simulations and it develops from the creeping flow 
solution as a function of the Reynolds number. For this case, in relation to a 
drag enhancement of ten percent the degree of asymmetry of the streamwise 
profile is of the same order. 

In figureini the relative change 1 — (u) / (u) is considered for C, = 0.25, C, = 0.5 
and C, = 0.75. The result is that using both LBM and FDM, excellent agree- 
ment is observed with respect to the functional behavior ~ i?e^ down to 
relative changes of 10~^. This is a stringent benchmark test for both numer- 
ical procedures for their degree of accuracy. 

Finally, in figure^! the value of F, as extracted from the numerical simula- 
tions is shown as a function of ^. As ^ approaches 1 the value of F is higher 
due to the fact that the effect of the non linear terms is enhanced and non 
trivially triggered by the boundary condition. 

6 Conclusions 

The importance of boundary slip velocity at finite Reynolds numbers has 
been studied for open channel flows. The local and global properties of an 
incompressible Newtonian fluid with mixed boundary conditions of no-slip 
and free shear have been characterized as a function of the boundary geome- 
try and Reynolds number. Although it would be easy to explain qualitatively 
the drag modification in the flow, its sensitivity with respect to the boundary 
condition and bulk properties may result a complicate subject. In our case, 
we have carried out numerical simulations with Lattice Boltzmann models 
(LBM) and Finite Difference methods (FDM) to highlight these properties 
and a main role is played by non linear effects that are induced by the non 
homogeneity of the boundary condition. For small Reynolds numbers we 
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have found interesting flow behaviors triggered by the boundary condition 
and very well controlled by the bulk Reynolds number. The overall physical 
picture is that a drag modification is induced by the non-uniform profile of 
the slip velocity at the boundary and it develops, for a given boundary real- 
ization, with a scahng law behavior as Re'^. The corresponding flow field is 
still laminar but looses its symmetry properties in the channel, with macro- 
scopic effects up to the ten percent in the effective slip length. The interplay 
and the physical correctness of both numerical schemes is then validated in 
the framework of the theoretical predictions for slip conditions realized on 
strips and reveals a pronounced efficiency of LBM with respect to FDM. This 
is of particular interest especially when we move towards the study of the 
same kind of problems in 3d geometries with much more realistic boundary 
conditions of slip spots and slip patterns on heterogeneous planar surfaces. 
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7 Appendix 

In this appendix we will detail the calculation leading to the analytical ex- 
pressions used throughout the text. Let us assume < i?e^ <^ 1 and expand 
u and p with respect to Rcr in the following form 

u = Relu^^^ + Relu^^^ + Re^u^"^^ + (40) 

p = Relp^^^ + i?ey^) + i?ey2) + ...^ (41) 

where the overscript (0), indicate the creeping flow solution, 

u^ReM^\ p^p^°\ (42) 

All the even orders (i.e., 0(i?e^")) of the velocity vector are identically zero 
because there is not any advection or driving forcing to balance with the 
divergence of the stress tensor containing the 0{Re'^) velocity. 
The continuity equation is now expressed as 

V-w^'^^^O (n = 0,l,...). (43) 
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The 0[Re^), 0{Re^) and 0(i?e^) momentum equations (fT^ are respectively 
written as 

= - Vp(°) + V^ii^o) + e,, 
V-(n(°)ii(°)) = -Vp« + VVi), 

For n = 0, 1 and 2, the boundary conditions are given by 

n(") = on NSW, 

(n-5("))xn = on FSW, 

where S*^"-* = | jVii*^") + (Vn^"))'^}. With these expansions, starting from 
equation (j25j) . we obtain the following relation 

{u)-{u) =Rel{u^^^) + Rel{u^^^) 

(44) 

The first and second terms in the rhs of equation ()44j] are zero because the 
kernel of the volume integral is written in the divergence form which can be 
rewritten using the area integral with no wall-normal flux: 




The third term in the rhs of (j44p is the only one that can take can take a 
non-zero value. Therefore, we obtain 

{u) - {u) Rel (n(i) ■ • V)u^'^^y]) as Re^ ^ +0. 
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It should be noted that for the expression 



{u) - {u) = TRel, 

the prefactor 

r = (n«-{(w(o)-v)w(°)}> 

is dependent upon the perturbed velocity component but independent of the 
Reynolds number. 
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Figure 1: The computational scheme for the numerical simulations. We use 
a 2d channel with the same boundary conditions on the top and lower walls. 
Free shear conditions are imposed on a width H on both top and lower 
walls giving a slip percentage E, = H/L^.. We also impose periodic boundary 
conditions in the streamwise (x) direction. 
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Time steps [xlO"^] 

Figure 2: Time evolution of the total kinetic energy, {ul+Uy), as a function of 
time steps. In the main figure we show with friction Reynolds number 

Rcr = 2.245 and parameters = U,Ly = 84,i/ = 1/6 -dP/dx = 1.88910"^ 
and ^ = 0.5. The agreement between LBM (o) and FDM (straight lines) is 
very good and at the same time the LBM converges more rapidly (~ 5 times 
faster). Inset: we report the same kind of plot for the case Rcr = 10.04, 
= 64,Ly = 168,z/ = 1/Q-dP/dx = 4.7210"^ and ^ = 0.5. Again LBM 
and FDM agree very well and in this case LBM is up to 15 times faster than 
FDM. 
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Figure 3: Velocity profiles along a channel with mixed boundary conditions 
of free shear and no-slip. In the inlet region (o for LBM and straight lines for 
FDM), due to the no slip condition the profile is almost a parabolic one with 
zero velocity at the boundaries and constant concavity, as predicted from the 
Poiseuille solution. In the middle of the channel (x for LBM and straight 
lines for FDM), the free shear condition is such to produce an enhancement 
of the local slip properties at the boundaries with a reduction of the profiles 
in the bulk region. In both cases all the profiles are normalized with respect 
to the center channel velocity in the Poiseuille limit. 




23 



0.30 



0.25 ■ 



0.20 h-J 



0.15 ■ 



0.10 ■ 



0.05 ■ 



0.00 








0.2 



0.4 



0.6 



0.8 



Figure 4: Slip length normalized to the pattern dimension as a function 
of the slip percentage ^. In the main figure we show a first order extrap- 
olation for LBM (o) and FDM (x). Both numerical schemes are com- 
pared with the analytical estimate (straight line) given in pH]- The nu- 
merical simulations are carried out with the following set of parameters 
Rcr = 2.245,L^ = QA,Ly = 84,z/ = 1/6 -dP/dx = l.SSgiQ-^.Inset: the 
same as the main figure but for the case of a fourth order approximation. In 
Both figures, being the Reynolds number small but finite, small discrepancies 
are observed between the numerics and the analytical results. 
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Figure 5: Drag modification at finite Reynolds numbers. We show the relative 
departure from the creeping flow solution, 1 — {u)/{u) (o), as a function of 
the friction Reynolds number Rcr- To verify the correctness of relation (jH^ 
we also plot the right hand side of this relation (x). All the results have 
been obtained with FDM with numerical grid mesh = Q4:,Ly = 64 and 
parameters u = 1/6,^ = 0.5. The pressure gradient dP/dx has then been 
changed to vary i?e^. To emphasize the scaling behavior with respect to -Re,-, 
the power law Re^ is also plotted. 
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Figure 6: Slip length at finite Reynolds numbers. We show the volume av- 
eraged slip length normalized to the creeping flow counterpart, as a function 
of the friction Reynolds number Rct-. All the results have been obtained 
with FDM with numerical grid mesh = 64, Ly = 64 and parameters 
u = 1/6,^ = 0.5. The pressure gradient dP/dx has then been changed to 
vary Rcr- In order to highlight the scaling behavior with respect to Rcr-, the 
power law behavior i?e^ is also represented. Note that for ~ 10, the 
overall slip length differs from its creeping flow counterpart of the order of 
ten percent. 
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Figure 7: Streamwise profiles for small Reynolds number. We plot the 
streamwise profiles as a function of the relative position along the channel 
for different distances from the wall (?/„_,_i with n = 0, 1,2, ... and the wall 
being located at yo = 0). Here we show with friction Reynolds number 

Rcr = 2.245 and parameters = QA,Ly = 84,i/ = 1/6 -dP/dx = 1.88910"^ 
and ^ = 0.5. The velocity is normalized with respect to the center channel 
velocity of the corresponding Poiseuille profile and both LBM (o) and FDM 
(straight lines) indicate an almost symmetric configuration with respect to 
the center of the slip area (see dotted line). This is predicted by the sym- 
metry properties of the Stokes solution {Rcr = 0) and it is here expected to 
hold with some very small corrections due to the finite Reynolds effects. 
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Figure 8: Same as figure[7|with the following set of parameters: RCr = 10.04 
and = QA,Ly = 168,i/ = 1/Q-dP/dx = 4.7210-^^ = 0.5. Now, due to the 
increase of the Reynolds number the profile becomes more asymmetric with 
respect to the case of figure To highlight this effect we plot the symmetry 
axis of the free shear strip to be compared with the symmetry properties of 
the flow. 
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Figure 9: Drag modification and its scaling behavior with respect to the 
friction Reynolds number for different values of the slip percentage ^. We 
show the relative departure from the creeping flow solution, 1 — {u) / {u) 
estimated with LBM (circles for C, = 0.25, squares for ^ = 0.5 and diamonds 
for ^ = 0.75) and FDM (straight lines). All the results have been obtained 
with FDM with numerical grid mesh = 64, Ly = 64 and parameters 
u = 1/6,^ = 0.5. The pressure gradient dP/dx has then been changed to vary 
Rct- In the small Re^ regime it is observed the scaling relation proportional 
to -Re^ predicted by our analytical approach. 
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